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ABSTRACT 

Accretion-powered X-ray pulsars are among the most luminous X-ray sources 
in the Galaxy. However, despite decades of theoretical and observational work 
since their discovery, no satisfactory model for the formation of the observed 
X-ray spectra has emerged. In particular, the previously available theories are 
unable to reproduce the power-law variation observed at high energies in many 
sources. In this paper, we present the first self-consistent calculation of the spec- 
trum emerging from a pulsar accretion column that includes an explicit treatment 
of the energization occurring in the shock. Using a rigorous eigenfunction expan- 
sion method based on the exact dynamical solution for the velocity profile in the 
column, we obtain a closed-form expression for the Green's function describing 
the upscattering of radiation injected into the column from a monochromatic 
source located at the top of the thermal mound, near the base of the flow. The 
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Green's function is convolved with a Planck distribution to calculate the radi- 
ation spectrum resulting from the reprocessing of blackbody photons emitted 
by the thermal mound. We demonstrate that the energization of the photons 
in the shock naturally produces an X-ray spectrum with a power-law shape at 
high energies and a blackbody shape at low energies, in agreement with many 
observations of accreting X-ray pulsars. 

Subject headings: methods: analytical — pulsars: general — radiation mecha- 
nisms: non-thermal — shock waves — stars: neutron — X-rays: stars 

1. INTRODUCTION 

Since the discovery of the first known pulsating X-ray sources Her X-l and Cen X- 
3 more than three decades ago (Giacconi et al. 1971; Tananbaum et al. 1972), over 50 
new sources have been detected in the galaxy and the Magellanic Clouds, with luminosities 
in the range L x ~ io 34-38 ergs s _1 and pulsation periods 0.1 s < P < 10 3 s. X-ray pulsars 
include a variety of objects powered by rotation or accretion, as well as several anomalous 
X-ray pulsars whose energy source is currently unclear. The emission from X-ray pulsars 
in binary systems is fueled by the accretion of material from the "normal" companion onto 
the neutron star, with the flow channeled onto one or both of the magnetic poles by the 
strong field. During the accretion process, gravitational potential energy is converted into 
kinetic energy, which escapes from the column in the form of X-rays as the gas decelerates 
through a radiative shock before settling onto the stellar surface. In the accretion-powered 
sources, which are the focus of this paper, the X-ray spectra are often well fitted using a 
combination of a power-law spectrum plus a blackbody component with a temperature in the 
range T ~ 10 6 ~ 7 K (e.g., Coburn et al. 2002; di Salvo et al. 1998; White et al. 1983). Most 
spectra also display quasi-exponential cutoffs at E ~ 20 — 30 keV, and there are indications of 
cyclotron features and iron emission lines in a number of sources. The observations suggest 
typical magnetic field strengths of ~ 10 12-13 G. 

Although accretion-powered X-ray pulsars are among the most luminous sources in the 
galaxy, previous attempts to calculate their spectra based on static or dynamic theoretical 
models have generally yielded results that do not agree very well with the observed profiles 
(e.g., Meszaros & Nagel 1985a,b; Nagel 1981; Yahel 1980; Klein et al. 1996). Hence there 
is still no clear understanding of the basic spectral formation mechanism in X-ray pulsars 
(see the discussion in Coburn et al. 2002). Given the lack of a viable theoretical model, 
X-ray pulsar spectral data have traditionally been fit using multicomponent forms that 
include absorbed power-laws, cyclotron features, iron emission lines, blackbody components, 
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and high-energy exponential cutoffs. The resulting parameters are sometimes difficult to 
relate to the physical properties of the source. Motivated by the lack of a comprehensive 
theoretical model for X-ray pulsar spectral formation, we reconsider here the physical picture 
originally proposed by Davidson (1973), in which the accreting gas passes through a radiative, 
radiation-dominated shock before settling onto the surface of the star. We illustrate the 
accretion/emission geometry schematically in Figure 1. Most of the photons emitted from 
the accretion column are produced in the dense "thermal mound" located inside the column, 
just above the stellar surface. The blackbody photons created in the mound are upscattered 
in the shock, and eventually diffuse through the walls of the column. The escaping photons 
carry away the kinetic energy of the gas, thereby allowing the plasma to settle onto the 
surface of the star. 



1.1. Bulk Comptonization 

The strong compression that occurs as the plasma crosses the radiative shock renders 
it an ideal site for first-order Fermi energization (i.e., "bulk" or "dynamical" Comptoniza- 
tion) of the photons produced by the thermal mound. In the bulk Comptonization process, 
particles experience a mean energy gain if the scattering centers they collide with are in- 
volved in a converging flow (e.g., Laurent & Titarchuk 1999; Turolla, Zane, & Titarchuk 
2002). By contrast, in the thermal Comptonizaton process, particles gain energy due to 
the stochastic motions of the scattering centers via the second-order Fermi mechanism (e.g., 
Sunyaev & Titarchuk 1980; Becker 2003). In the X-ray pulsar application, the scattering 
centers are infalling electrons, and the energized "particles" are photons. Since the inflow 
speed of the electrons in an X-ray pulsar accretion column is much larger than their thermal 
velocity, bulk Comptonization dominates over the stochastic process except at the highest 
photon energies (Titarchuk, Mastichiadis, & Kylafis 1996), as discussed in § 7. The failure of 
the current models to generate the power-law spectral shape characteristic of X-ray pulsars 
probably stems from the neglect of the critical contribution of the shock in upscattering the 
soft radiation produced in the thermal mound (e.g., Burnard, Arons, & Klein 1991). We 
demonstrate in this paper that bulk Comptonization of the radiation in the accretion shock 
can produce spectra very similar to the power-law continuum seen in many X-ray pulsars. 
The main results presented here were summarized by Becker & Wolff (2005). 

Dynamical Comptonization has been considered by a number of previous authors in a va- 
riety of astrophysical situations, including spherical inflows (Payne & Blandford 1981; Colpi 
1988; Schneider & Bogdan 1989; Titarchuk, Mastichiadis, & Kylafis 1997), accretion disks 
around black holes (Laurent & Titarchuk 2001; Titarchuk, Cui, & Wood 2002; Titarchuk & 
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Shrader 2002), spherical outflows (Becker & Begelman 1986; Titarchuk, Kazanas, & Becker 
2003), and the transport of cosmic rays in the solar wind (Parker 1965; Stawicki, Fichtner, 
& Schlickeiser 2000). The process was also studied in the context of neutron star accretion 
by Titarchuk, Mastichiadis, & Kylafis (1996) and Mastichiadis & Kylafis (1992), although 
their results are not applicable to X-ray pulsar accretion columns due to the assumption of 
spherical symmetry and the imposition of a power-law velocity profile. Lyubarskii & Sunyaev 
(1982) analyzed dynamical Comptonization in the context of plane-parallel pulsar shocks, 
but the velocity profile they utilized is not consistent with the dynamics of X-ray pulsar ac- 
cretion, and furthermore they did not account for the escape of radiation through the walls 
of the column. The relevance of the bulk Comptonization process for spectral formation in 
X-ray pulsars was also recognized by Burnard, Arons, & Klein (1991), who pointed out that 
thermal spectra alone are much too soft to explain the observed emission, and suggested 
that upscattering in the shock could provide the required hardening. 

The fundamental character of the Green's function describing both thermal and bulk 
Comptonization was studied by Titarchuk & Zannias (1998) for the case of accretion onto a 
black hole. These authors established that the Green's function can be approximated using 
a broken power-law form with a central peak between the high- and low-energy portions of 
the spectrum, for either type of Comptonization. Furthermore, they concluded that bulk 
Comptonization dominates over the thermal process if the electron temperature T < 10 7 K. 
The theory developed in the present paper represents an extension of the same idea to 
neutron star accretion, resulting in the exact Green's function for the X-ray pulsar spectral 
formation process. In agreement with Titarchuk & Zannias (1998), we find that the bulk 
Comptonization process is dominant for the temperature range relevant for X-ray pulsars. 
However, the presence of the event horizon in the black-hole application treated by these 
authors makes their problem fundamentally different from the neutron star accretion problem 
studied here, because the neutron star obviously possesses a solid surface. Our work therefore 
represents the first exact, quantitative analysis of the role of bulk Comptonization in the X- 
ray pulsar spectral formation process. 



Radiation pressure governs the dynamical structure of the accretion flows in bright 
pulsars when the X-ray luminosity satisfies (Becker 1998; Basko & Sunyaev 1976) 



1.2. 



Radiation-Dominated Flow 




ergs s 



-i 




- 5 - 



where r is the polar cap radius, M* and R* denote the stellar mass and radius, respectively, 
<j t is the Thomson cross section, and cry and a±_ represent the electron scattering cross 
sections for photons propagating parallel or perpendicular to the magnetic field, respectively. 
When the luminosity of the system is comparable to L crit , the radiation flux in the column 
is super-Eddington and therefore the radiation pressure greatly exceeds the gas pressure 
(Becker 1998). In this situation the gas passes through a radiation-dominated shock on 
its way to the stellar surface, and the kinetic energy of the gas is carried away by the high- 
energy radiation that escapes from the column. The strong gradient of the radiation pressure 
decelerates the material to rest at the surface of the star. The observation of many X-ray 
pulsars with L x ~ 10 36 " 38 ergs s" 1 implies the presence of radiation-dominated shocks close 
to the stellar surfaces in these systems (White et al. 1983; White et al. 1995). Note that 
radiation-dominated shocks are continuous velocity transitions, with an overall thickness of 
a few Thomson scattering lengths, unlike traditional (discontinuous) gas-dominated shocks 
(Blandford & Payne 1981b). 

In luminous X-ray pulsars, the pressure of the radiation determines the dynamics of the 
flow, while the dynamics in turn determines the shape of the radiation spectrum. Hence 
the flow dynamics and the radiative transport are intimately connected, which makes this 
a complex and nonlinear "photohydrodynamical" problem. One implication is that the 
photons scattering through the column cannot be regarded as "test particles" since it is their 
pressure that dictates the structure of the flow. We must therefore solve for the radiation 
spectrum and the velocity profile in a self-consistent manner. Radiation pressure may also 
play an important role in the dynamics of moderate-luminosity pulsars due to the strong 
dependence of the electron scattering (cyclotron) cross section on the magnetic field strength 
(Langer & Rappaport 1982). 

Becker (1998) and Basko & Sunyaev (1976) have considered the dynamics of radiation- 
dominated pulsar accretion flows, including the effect of the shock and the escape of radiation 
energy from the column. They find that the gas is decelerated to zero velocity in a manner 
consistent with accretion onto a solid body, but they do not consider the shape of the es- 
caping radiation spectrum. Conversely, the dynamics and the radiation spectrum have been 
worked out self-consistently by Blandford & Payne (1981b) for the case of a "standard" 
(Rankine-Hugoniot) radiation-dominated shock, which does not include radiation escape. In 
the Blandford & Payne approach, the radiation spectrum everywhere in the flow is deter- 
mined by solving the transport equation using the incident radiation field as a boundary 
condition. The solutions have a power-law character at high energies, with a spectral index 
that depends on the Mach number of the upstream flow. Since the Rankine-Hugoniot shocks 
studied by Blandford & Payne have a conserved energy flux, the downstream velocity is never 
less than one-seventh its upstream value. Consequently the associated spectral solutions are 
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Fig. 1. — Schematic depiction of gas accreting onto one of the magnetic polar caps of a 
neutron star. 

not applicable to the X-ray pulsar case, since pulsar shocks must be radiative in nature so 
that the kinetic energy of the infalling material can be removed and the flow brought to rest 
at the stellar surface. 

The analytical approach employed here parallels the work of Blandford & Payne (1981b), 
except that the velocity profile they utilized is replaced with the appropriate solution for 
an X-ray pulsar accretion column, and the loss of radiation energy from the column is 
incorporated into the transport equation using an escape-probability formalism. In the 
spirit of the Blandford & Payne study of spectral formation in radiation dominated Rankine- 
Hugoniot shocks, we are interested here in exploring the direct effect of the radiative accretion 
shock on the spectrum emerging from an X-ray pulsar. Hence it is not our goal in this paper 
to develop complete models that include additional processes such as free-free emission and 
absorption, cyclotron features, and iron emission lines. However, even without including 
any of these processes, we are able to demonstrate qualitative agreement with X-ray pulsar 
spectra. This suggests that energization in the accretion shock is one of the most important 
aspects of spectral formation in X-ray pulsars. 

The remainder of the paper is organized as follows. In § 2 we briefly review the conser- 
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vation equations and discuss the resulting dynamical solution for the velocity profile of the 
accreting gas. The transport equation governing the scattering of the radiation inside the 
accretion column is analyzed in § 3, and various constraints are derived in order to ensure 
that the radiative transfer is modeled in a manner consistent with the flow dynamics. In 
§ 4 we obtain the exact analytical solution for the Green's function and convolve it with 
the Planck distribution emitted at the thermal mound to compute the radiation spectrum 
emerging from the accretion column. In § 5 we present results obtained for the count-rate 
spectra using parameters corresponding to specific pulsars, and we compare the theoretical 
spectra with the observational data. The implications of our results for the production of 
X-ray spectra in pulsars are discussed in § 7. Additional mathematical details are provided 
in a series of Appendices. 



2. DYNAMICS OF RADIATION-DOMINATED FLOW 

The model analyzed here is based on the dynamical picture suggested by Davidson 
(1973), which is illustrated schematically in Figure 1. The accretion scenario corresponds 
physically to the flow of a mixture of gas and radiation inside a magnetic "pipe" that is 
sealed with respect to the gas, but is transparent with respect to the radiation. The accretion 
column incorporates a radiation-dominated, radiative shock located above the stellar surface. 
Although the shock is extended and the velocity profile is continuous, the flow possesses a 
well-defined sonic surface, as indicated in Figure 1. The model also includes a dense thermal 
mound located at the base of the flow, where local thermodynamic equilibrium prevails. Due 
to its blackbody nature, the thermal mound acts as both a photon source (emitting a Planck 
spectrum) and a photon sink (absorbing all incident radiation). Hence the surface of the 
mound represents the "photosphere" for photon creation and absorption, and the opacity is 
dominated by electron scattering above this point. 

In our approach to the problem, we assume that the upstream flow is composed of 
pure, fully ionized hydrogen gas moving at a highly supersonic speed, which is the standard 
scenario for accretion-powered X-ray pulsars (Basko & Sunyaev 1975, 1976). Our transport 
model employs a cylindrical, plane-parallel geometry, and therefore the velocity, density, 
and pressure are functions of the distance above the stellar surface, but they are all constant 
across the column at a given height. In this situation, the vertical variation of the flow 
velocity is given by the exact analytical solution obtained by Becker (1998) and Basko & 
Sunyaev (1976), which describes the settling of the transonic flow onto the surface of the star. 
A detailed consideration of the angular and energy dependences of the electron scattering 
cross section is beyond the scope of the present paper. We shall therefore follow Wang 
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& Frank (1981) and Becker (1998) by treating the directional dependence of the electron 
scattering in an approximate way in terms of the energy-averaged cross sections a\\ and a±, 
describing respectively the scattering of photons propagating either parallel or perpendicular 
to the magnetic field. The gas in the accretion column is radiation-dominated, and therefore 
it is the pressure of the photons that decelerates the infalling plasma to rest. We will use this 
fact in § 6.1 to evaluate the self-consistency of the coupled radiative transport/dynamical 
model based on the analytical solution for the radiation pressure profile obtained below. 



2.1. Conservation Equations 

The one-dimensional, time-dependent Euler equations governing the flow of a radiation- 
dominated gas with mass density p, flow velocity v, and photon energy density U inside a 
cylindrical accretion column are 

?E - -^L o\ 

dt~ dx' [ ] 

9 , N BI 

Ft = "S ' (3) 

— f - pv 2 + U j = -— + U csc + f/ abs + U emit , (4) 
where the mass, momentum, and energy fluxes are given respectively by 

J = pv , (5) 

I = P + pv \ ( 6 ) 

1 BP 

E = -pv 3 + Pv + Uv-c— , (7) 

2 &q 

with P = U/3 denoting the radiation pressure. The spatial coordinate x increases along the 
column axis, in the direction of the flow, and t\\ denotes the associated electron-scattering 
optical depth, which is related to x via 

dr\\ = n e cry dx , (8) 

where n e = p/m p is the electron number density, m p is the proton mass, and a\\ is the 
electron scattering cross section parallel to the magnetic field (i.e., the column axis). The 
two coordinates x and t\\ are calibrated so that they each vanish at the sonic point, as 
explained below. 

The quantities C/ esc , U&bs, and £7 C mit appearing in equation (4) represent the rates per 
unit volume at which radiation energy escapes through the column walls, or is absorbed or 
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emitted by the gas, respectively. When the system is radiation-dominated as assumed here, 
the gas itself has little thermal energy to contribute to the radiation field, although it does 
possess a large reservoir of bulk kinetic energy which is transferred to the photons through 
the first-order Fermi energization process in the shock. In this situation, the emission and 
absorption processes must nearly balance throughout the flow, and therefore we can write 

^abs + ^emit = , (9) 

in which case equation (4) reduces in a steady state to 

& = °- ■ (10) 

Although emission and absorption do not effect the flow dynamics in this situation, these 
processes profoundly influence the shape of the radiation spectrum, as discussed in §§ 3 and 4. 

Within the context of our one-dimensional picture, the rate at which radiation energy 
diffuses through the walls of the column can be approximated using an escape-probability 
formalism. If the column cross section is optically thick to electron scattering, as expected, 
then the rate at which radiation energy escapes through the walls is given by 

U , x 

C/esc = -— , (11) 
''esc 

where the mean escape time, t eS c, is given by 

t esc = — , w± = — , r± = n e a±r , (12) 

W± T_|_ 

with w± denoting the diffusion velocity perpendicular to the rr-axis, and t± representing the 
perpendicular optical thickness of the cylindrical accretion column with radius tq. Note that 
t± and tesc are functions of x through their dependence on n e . Becker (1998) confirmed that 
the diffusion approximation employed in equation (11) is valid because r± > 1 for typical 
X-ray pulsar parameters. Equations (5) and (12) can be combined to express the mean 
escape time as 

tesc = • (13) 

rripcv 

Since the escape timescale is inversely proportional to the flow velocity, the column becomes 
completely opaque at the surface of the neutron star due to the divergence of the electron 
number density there. The relationship between the escape-probability approximation and 
the physical distribution of radiation inside the accretion column is further discussed in § 6.3. 
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2.2. Dynamical Solution 



In the steady-state situation of interest here, the mass flux J and the momentum flux I 
are both conserved, but the energy flux E decreases as the gas accretes onto the surface of 
the neutron star due to the emission of radiation through the column walls. The momentum, 
energy, and mass conservation equations can be combined in this case to show that the flow 
velocity v satisfies the second-order nonlinear differential equation (Becker 1998) 



4/,) , (14) 



where 



v c v 

T = T|| — , fi= — 



(15) 

c v c 

Here, v c represents the flow velocity at the sonic point, where we set x = r = and M. — 1, 
with the radiation Mach number Ai defined by 

v 2 _ 4 P 

c s 3 p 



M 



(16) 



The quantity c s denotes the speed of sound in the radiation-dominated gas. In the far 
upstream region, we assume that M. — > oo which is an excellent approximation in pulsar 
accretion flows (Basko & Sunyaev 1975, 1976). 

Becker (1998) showed that in order for the flow to come to rest at the stellar surface, 
the parameters r , J, cr±, and an must satisfy the relation 



nip c 2 



4 
3 



(17) 



The corresponding exact solution for the velocity v as a function of r in the steady-state 
situation of interest here is given by 



v(t) 



7v, 



1 — tanh 



where 



2 i i 1 
- tanh" 1 - 

7 V 7 



(r - n) 



0.041 . 



(18) 



(19) 



Note that fi = 1 at the sonic point (r = 0) as required. The velocity can also be written as 
an explicit function of the position x using (Basko & Sunyaev 1976) 



v(x) 



7v r 



-l+x/Xs 



(20) 
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where x st is the distance between the sonic point and the stellar surface, which can be 
evaluated using equation (4.16) from Becker (1998) to obtain 



** = iTfUJ ln UJ' (21) 

Note that x increases in the direction of the flow, which comes to rest at the surface of the 
star as expected. The height above the stellar surface h is related to the coordinate x via 

h(x) = Xgt — x ■ (22) 

The incident (upstream) velocity of the infalling material is expected to be close to the 
free-fall velocity onto the stellar surface, v s , given by 

/2GMA" 2 

%=(— j ■ (23) 

Since v — * (7/4) v c in the upstream region according to equation (20), it follows that v c = 
(4/7) Vtf , and therefore the flow velocity at the sonic point is given by 

V < = -7{—) ■ <24) 

With the velocity solution given by equation (20), the corresponding pressure profile 
can be obtained by combining equations (5) and (6), which yields 

P(x) = I - Jv(x) . (25) 

The pressure increases to a maximum value at the surface of the star, where v — > 0. Based 
on equations (5) and (16), we can reexpress the Mach number as 

Since the Mach number approaches infinity in the far upstream region and v — > (7/4) v c , it 
follows that P must vanish asymptotically and consequently the upstream momentum flux 
is dominated by the ram pressure of the freely falling matter. We therefore obtain 

I = Jv s . (27) 

By combining equations (24), (25) and (27), we find that the pressure profile is related to 
the velocity variation by 

7 v(x) 



P(X) = Jv r 



(28) 
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Using equation (20) to eliminate the velocity v(x) yields the closed-form result 

P{x) = Lj Vc ^ . (29) 

The pressure vanishes in the upstream limit (x — > — oo) as expected, and at the surface of 
the star (x = x st ), the pressure achieves the stagnation value 

P st = P(x st ) = ^Jv c . (30) 

3. RADIATIVE TRANSFER IN THE ACCRETION COLUMN 

If the gas is radiation-dominated and fully ionized, then the photons interact with 
the matter primarily via electron scattering, which controls both the spatial transport and 
the energization of the radiation. In this situation the opacity is dominated by electron 
scattering, and therefore absorption is negligible, except at the surface of the thermal mound 
located inside the cylindrical accretion column at x = xq. We are interested in obtaining 
the photon distribution, f(xo,x,e), measured at position x and energy e resulting from 
the reprocessing of blackbody radiation emitted from the mound. The normalization of / 
is defined so that e 2 f(x , x, e) de gives the number density of photons in the energy range 
between e and e + de, and therefore f = 8nn/ (c 3 /i 3 ), where n is the occupation number. The 
self-consistency of the solution obtained for the radiation distribution / will be confirmed 
by verifying that the associated pressure distribution agrees with the dynamical result given 
by equation (29). 

3.1. Transport Equation 

In the cylindrical geometry employed here, the photon distribution / satisfies the trans- 
port equation (e.g., Blandford & Payne 1981a; Becker 1992; Parker 1965; Becker & Wolff 
2005) 

df _ v df | dvedf | d t c df\ | S(e) f 
dt dx dx 3 de dx \3n e <7|| dx J -nr\ t csc 

where e is the photon energy, x is the spatial coordinate along the column axis, and v = v(xq) 
is the flow velocity at the top of the thermal mound. The terms on the right-hand side of 
equation (31) represent advection, first-order Fermi energization ("bulk Comptonization" ) 
in the converging flow, spatial diffusion parallel to the column axis, the blackbody source, 
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escape of radiation from the column, and the absorption of radiation at the thermal mound, 
respectively. We are interested here in the steady-state version of equation (31) with df/dt = 
0. The photon number and energy densities associated with the distribution / are given 
respectively by 



In the present application, we are primarily interested in exploring the effect of the dynamics 
itself (i.e., the radiative shock) on the observed X-ray emission. The transport equation (31) 
therefore does not include additional processes such as thermal Comptonization or cyclotron 
emission and absorption that are likely to be important in X-ray pulsars. In general, the 
neglect of thermal Comptonization is reasonable because at the temperatures typical of X- 
ray pulsars (T ~ 10 6 ~ 7 K), the kinetic energy associated with the bulk flow far surpasses the 
thermal energy. However, the failure to include the electron recoil associated with thermal 
Comptonization renders the model unable to reproduce the quasi-exponential cutoffs often 
observed in X-ray pulsar spectra at high energies. Furthermore, in some sources thermal 
Comptonization appears to be necessary in order to flatten the spectrum in the energy 
range e ~ 5 — 20 keV. We provide further discussion of this issue in § 7. 

When integrated over the photon energy spectrum, the first-order Fermi process con- 
sidered here corresponds to the PdV work done on the radiation by the compression of 
the background plasma as it accretes onto the stellar surface. The mean escape time, t esc , 
is given by equation (13), and the source term S(e) is defined so that e 2 S(e) de represents 
the number of photons injected into the accretion column per second in the energy range 
between e and e + de from the blackbody surface at x = x$. The absorption term in equa- 
tion (31) containing the dimensionless constant f3 must be included due to the blackbody 
nature of the thermal mound, which acts as both a source and a sink of radiation. In the 
radiation-dominated situation considered here, essentially all of the pressure is provided by 
the photons, and therefore the rate of absorption of radiation energy at the top of the mound 
must equal the energy emission rate when integrated with respect to the photon frequency 
(see eq. [9]). We will use this energy balance argument to compute a self-consistent value 
for (3 below. 

Next we need to specify the appropriate form for the source term S(e) in order to treat 
the injection of the blackbody radiation from the surface of the thermal mound. The flux per 
unit energy measured at the surface of an isotropically radiating object is equal to ix times 
the intensity (Rybicki & Lightman 1979). It therefore follows that the amount of energy 
emitted per second from the upper surface of the thermal mound (with area ttTq) in the 




(32) 



-14- 



energy range between e and e + de is given by 

e 3 S(e) de = ixr\ • n B e (e) de , (33) 

where 

2e 3 1 

B e (e) = —j (34) 

v ' h 3 c 2 e £ / fcT ° - 1 

denotes the blackbody intensity and To = T(xq) is the gas temperature at the surface of the 
mound. Note that the units for B e are ergs s _1 ster -1 cm -2 erg -1 . Hence we obtain for the 
source term 

2n 2 r 2 1 

S{e) = ~hF^ e'/wo - 1 ' ( 35 ) 
and consequently the total energy injection rate is given by 

POO 

/ e 3 S (e) de = n a Tq oc ergs s _1 , (36) 
Jo 

where a is the Stephan-Boltzmann constant. 



3.2. Energy Balance at the Thermal Mound 

In a radiation-dominated X-ray pulsar accretion column, the emission of fresh blackbody 
radiation energy at the surface of the thermal mound is almost perfectly balanced by the 
absorption of energy, as expressed by equation (9). Most of the energy appearing in the 
emergent X-rays is therefore not supplied from the internal energy of the gas, but rather 
from its bulk kinetic energy, which is transferred directly to the photons via collisions with 
infalling electrons. We can gain some insight into the energy balance at the thermal mound 
by rewriting the steady-state version of the transport equation (31) in the flux-conservation 
form (Skirling 1975; Gleeson & Axford 1967) 

dH e 1 d f 3 df\ S(e) f 

^ = -^e{ eV ^) + ^f 6{x - Xo) -IZ- (3Vo6{x - Xo)f > (37) 



where 



^-T^-TT-TTT < 38 > 

3 n P an ox 3 oe 



represents the "specific flux" (Becker 1992). By operating on equation (37) with f^°e 3 de 
and utilizing equations (32) and (36), we arrive at the photon energy equation 

~d^x = 1 ~d~x ~ T~ ~ ^ Wo( ^ x ~~ x °) U + aT oK x ~ x o) , (39) 
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where the integrated photon energy flux, Q, is given by 

f°° 4 c dU 

Q = / e 3 H e de = -vU — oc ergs cm' 2 s _1 . (40) 

J 3 3n e <T|| dx 

Equation (39) states that the divergence of the radiation energy flux Q is equal to the net 
rate of change of the photon energy density U due to the combined influence of compression, 
escape, absorption, and photon injection. Note that Q is related to the total energy flux E 
via (see eq. [7]) 

E = Q+ l -pv z . (41) 

Since the absorption and emission of radiation energy must balance at the surface of the 
mound in the radiation-dominated situation of interest here, we can integrate equation (39) 
with respect to x across the mound location to obtain 

lim Q(x + e) - Q(x - e) = -f3 v U + a T 4 = , (42) 

e^O 

where Uq = U(xq) denotes the radiation energy density at the mound surface. The fact that 
the photon energy flux Q is continuous at x = xq therefore requires that the dimensionless 
absorption constant (3 in the transport equation (31) be given by 

rrT 4 

P = ZTTT ■ (43) 

By utilizing this expression for /3, we ensure that the results obtained for the radiation spec- 
trum J'(xq, x, e) are consistent with the dynamical structure of the accretion column. This is 
confirmed after the fact by calculating the radiation pressure P — U/3 using equation (32) 
and comparing the result with the dynamical solution given by equation (29). 



4. EXACT SOLUTION FOR THE RADIATION SPECTRUM 

In our approach to solving equation (31) for the spectrum J'(xq, x, e) inside the accretion 
column, we shall first obtain the Green's function, f G (xo, x, eo, e), which is the radiation 
distribution at location x and energy e resulting from the injection of N photons per second 
with energy eo from a monochromatic source at location xq. The determination of the Green's 
function is a useful intermediate step in the process because it provides us with fundamental 
physical insight into the spectral redistribution process, and it also allows us to calculate the 
particular solution for the spectrum / associated with an arbitrary continuum source S(e) 
using the integral convolution (Becker 2003) 

f(x , x, e) = f /G(X0 ^ e0,£) e 2 S(e ) de . (44) 
Jo vv 
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Note that the upper bound of integration is set equal to e because in our model all of 
the injected photons gain energy via the first-order Fermi process, and none lose energy. 
The technical procedure used to solve for the Green's function involves the derivation of 
eigenfunctions and associated eigenvalues based on the set of spatial boundary conditions 
for the problem (see, e.g., Blandford & Payne 1981b; Payne & Blandford 1981; Schneider & 
Kirk 1987; Colpi 1988). 

The steady-state transport equation governing the Green's function is (cf. eq. [31]) 

M = * £ ?L + » (_£_ f.) + ^-^f-^ ^p Vo S{x . Xo) fa , (45) 
Ox dx 3 Oe ox \6n e a\\ ox J 7rr o e o ^esc 

and the associated radiation number and energy densities are given by (cf. eqs. [32]) 

poo pco 

n G (x)= I e 2 f G (x ,x,e ,e)de , U G (x) = / e 3 / G (x , x, e , e) de . (46) 
Jo Jo 

Further simplification of the mathematical derivation is possible if we work in terms of the 
new spatial variable y, defined by 

rj\ -l+x/x s t 



y{x) -^j • (47) 

where x st is given by equation (21). Note that y — > in the far upstream region (x — > — oo), 
and y — > 1 at the surface of the star (x — > x st ). 

Based on equations (20) and (47), we find that the variation of the velocity v as a 
function of the new variable y is given by the simple expression 

v(y) = 7 -^(l-y). (48) 

Likewise, we can also combine equations (29) and (47) to show that the exact dynamical 
solution for the radiation pressure profile as a function of y is given by 

P{y) = \jv c y. (49) 

Note that in the limit y — > 1, the pressure approaches the stagnation value P st = (7/4)Jv c 
in agreement with equation (30). Equations (29) and (49) will prove useful when we seek to 
confirm the validity of the results obtained for the photon distribution f(xo, x, e) in § 5. 

Utilizing equations (5), (13), (17), and (48) along with the differential relation 

dx r f(T±\ 1/2 _x , s 

5TwsUJ v ' (50) 
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we find that equation (45) can be transformed from x to y to obtain 



dy 2 \ 4 J dy 4 de \ 4y 

3(3v 5(y - y ) f G 3 N 5(e - e ) 5(y - y ) 



7v c lnrle\ v c 



(51) 



where we have made the definition 

y = y(z ) = ( - ) , (52) 



rj\ —1+X0/X s t 



so that y denotes the value of y at the top of the thermal mound (see eq. [47]). According 
to equation (48), the flow velocity at the top of the mound, vo, is related to yo via 

* = (53) 

Note that we can now write the Green's function as either f G (xo, x, eo, e) or f G (yo,y,eo,e) 
since (x,x ) and (y,yo) are interchangeable via equations (47) and (52). 



4.1. Separation Solutions 

When e > e , equation (51) is separable in energy and space using the functions 

f x (e,y) = e~ x g(\,y) , (54) 
where A is the separation constant, and the spatial function g satisfies the differential equation 



l-5y\dg fXy + y-1 



dy 



+ 



4y 



9 



3pv 5(y - y ) 



9 



(55) 



In order to avoid an infinite spatial diffusion flux &ty = yo, the function g must be continuous 
there, and consequently we obtain the condition 



lim g(X, y + e) - g(X, y - e) = 

£^0 



(56) 



We can also derive a jump condition for the derivative dg/dy at the top of the mound by 
integrating equation (55) with respect to y in a very small region around y = yo. The result 
obtained is 



lim — 

e^o dy 



y=yo+e 



dg_ 

dy 



y=yo-e 



3j3 
4 yo 



g(\yo) , 



(57) 
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where we have also used equation (53). The spatial eigenfunctions for our problem are 
denoted by 

g n (y) = g(K,y) , (58) 

where \ n represents the nth eigenvalue. These functions satisfy the continuity and derivative 
jump conditions given by equations (56) and (57), respectively, as well as the boundary 
conditions discussed below. 

The homogeneous version of equation (55) for g obtained when y ^ y has fundamental 
solutions given by 

<Pi(Ky) = y F (a, b; c; y) , (59) 
<pl(\,y) = y- 1/4 F(a - 5/4, b - 5/4; 2 - c; y) , (60) 

where F(a,b;c;z) denotes the hypergeometric function (Abramowitz & Stegun 1970), and 
the parameters a, b, and c are defined by 

9 - V17+ 16 A , 9 + V17 + 16 A 9 , , 

a = ^ , b= , c = - . 61 

8 8 ' 4 y J 

The "seed" photons injected into the flow from the thermal mound located near the base 
of the column are unable to diffuse very far up into the accreting gas due to the extremely 
high speed of the inflow. Most of the photons therefore escape through the walls of the 
column within a few scattering lengths of the mound, forming a "fan" type beam pattern, 
as expected for accretion-powered X-ray pulsars (e.g., Harding 1994, 2003). Based on this 
physical picture, we conclude that the eigenfunction g n must vanish in the upstream limit, 
y — > 0. Conversely, in the downstream limit, y — > 1, the gas settles onto the surface of the 
star and the radiation pressure achieves the stagnation value given by equation (30). Since 
this is a finite pressure, we conclude that g n must approach a constant as y — > 1. 

Analysis of the asymptotic behaviors of the hypergeometric functions <fi and ip\ shows 
that the first function vanishes and the second diverges as y — > 0. Hence in the region 
upstream from the thermal mound (y < yo), the eigenfunction g n must be given by ipi. 
However, in the downstream region (y > y ), the situation is more complicated because the 
two functions ipi and tp\ each diverge logarithmically as y — > 1. We must therefore utilize a 
suitable linear combination of these functions in order to obtain a convergent solution in the 
downstream region. Consequently we define the new function (see Appendix A for details) 

M s r ( c)r(?-t) M - r ( 2- J)r (a ) M • (62) 

which remains finite in the limit y — > 1 as required. 
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4.2. Eigenf unctions and Eigenvalues 

By utilizing the various relations derived above, we find that the global solution for the 
spatial eigenfunction g n {y) (eq. [58]) can now be written as 

9n{y) = n / \ x ^ (63) 

where the constant B n is determined using the continuity condition (eq. [56]), which yields 

B„ = ^4 . (64) 

We can combine equations (57), (63), and (64) to show that the eigenvalue equation for A is 
given by 

d(f 2 d(fi! 3/3</?i</?2 



<Pi -r; V?2 



= . (65) 

y=yo 



dy dy Ay 

The left-hand side of this expression can be evaluated using the exact solution for the Wron- 
skian (see Appendix B), 

W ^ y) = ^-^ = l V{a)V{2-c) —y ' (66) 

which is applicable for arbitrary values of A and y. Equations (65) and (66) can be combined 
to obtain the alternative form for the eigenvalue equation, 

[{l "' U] !h ' ] -/^i(A,y )^(A,y ), (67) 



3 Y{a) T(2 - c) 1 - y 

where a, b, and c are functions of A given by equations (61). The roots of this expression are 
the eigenvalues A = A n , and the associated eigenfunctions are evaluated using equation (63). 

The first eigenvalue, A , is especially important because it determines the slope of 
the high-energy spectrum emerging from the accretion column according to equation (54). 
The spectral index of the emitted photon count rate distribution, ao, is related to Ao via 
«o = Ao — 2. In Figure 2 we plot the photon index ao as a function of the dimensionless 
parameters (3 and yo- Note that «o is a double-valued function of yo for fixed /3, which is 
a consequence of the imposed velocity profile (eq. [48]). Physically, this behavior reflects 
the fact that it is always possible to achieve a desired amount of compression (first-order 
Fermi energization) by placing the source in a specific location in either the upstream or 
downstream region of the flow. We also observe that if we increase the absorption parameter 
(3 while holding y fixed, then a increases monotonically, and therefore the high-energy 
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Fig. 2. — High-energy power-law photon spectral index a plotted as a function of the source 
location y for the indicated values of the absorption parameter f3. Note the steepening of 
the radiation spectrum that occurs when (5 is increased for a fixed value of yo, which reflects 
the decreasing residence time for the photons in the plasma due to the enhanced absorption. 

spectrum becomes progressively steeper. This behavior is expected physically because as 
the absorption parameter is increased, the injected photons spend less time on average being 
energized by collisions with electrons before either escaping from the column or being ab- 
sorbed at the source location. The reduction in the amount of energization naturally leads 
to a steepening of the radiation spectrum. When (5 = 0, no absorption occurs, and the index 
a achieves it minimum (limiting) value of 2. This limit is unphysical, however, since it 
yields a divergent result for the total photon energy density U G according to equation (46). 



We demonstrate in Appendix C that the eigenfunctions g n (y) form an orthogonal set, as 
expected since this is a standard Sturm-Liouville problem. Once the eigenvalues and eigen- 
functions are determined, the solution for the Green's function can therefore be expressed 
as the infinite series 



4.3. Eigenfunction Expansion 




(68) 
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where the expansion coefficients C n are computed by employing the orthogonality of the 
eigenfunctions along with the condition 



f G (yo,y,e ,e) 



12 Nn 



e=e 



S(y-yo) , 



(69) 



which is obtained by integrating the transport equation (51) with respect to e in a very 
small range surrounding the injection energy eo- The result obtained for the nth expansion 
coefficient is 

i2 7V y " 3/4 ^(yo) 



7wrl v c I n 



where the quadratic normalization integrals I n are defined by 



I n = / y-^g n (y)dy 



(70) 



(71) 



In Appendix D, we show that the normalization integrals can be evaluated using the closed- 
form expression 

I n = K(X n ,y ), (72) 

where 

"*(a)+#(l-a) dlnipi d\mp 2 



and 



K(\,y)=3f3y- 3 /\l-y) V l(\,y) 



(17+ 16A)V2 
1 dT(z) 



(73) 



(74) 



T(z) dz 

This provides an extremely efficient alternative to numerical integration for the computation 
of I n . The eigenfunction expansion converges fairly rapidly and therefore we are able to 
obtain an accuracy of at least five significant figures in our calculations of f G by terminating 
the series in equation (68) after the first 20 terms. 



4.4. Green's Function for the Escaping Photon Spectrum 

Equation (68) represents the exact solution for the Green's function describing the 
radiation spectrum inside a pulsar accretion column resulting from the injection of iVo seed 
photons per unit time from a monochromatic source located at y = yo (or, equivalently, at 
x = x ). In the escape probability approach employed here, the associated Green's function 
for the photon spectrum emitted through the walls of the cylindrical column is defined by 

N^(x ,x,e ,e) = — - — f G (x , x, e , e) , (75) 

^esc 
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so that iV G dx de represents the number of photons escaping from the column per unit time 
between positions x and x + dx with energy between e and e + de. We remind the reader 
that the quantities (x,x ) and (y,yo) are interchangeable via equations (47) and (52), and 
therefore we are free to work in terms of the more convenient parameters (y, yo) without 
loss of generality. By substituting for t csc using equation (13) we can obtain the alternative 
expression 

• r n m p cv(y) e 2 

N e (yo,y,eo,e) = f G {y ,y,e ,e) . (76) 

J a j_ 

Eliminating J using equation (17) and substituting for v and f G using equations (48) and 
(68), respectively, yields the equivalent result 

7V e G (y , y, co, c) = (i - y) £ A, ( - ' ( 7? ) 

where the expansion coefficients D n are defined by 

r yo' I n \ a± / 

The Green's function iV G describes the photon distribution escaping from the accretion 
column as a function of energy e and location y for the case of monoenergetic photon in- 
jection. Analysis of iV G therefore reveals some interesting details about the energization of 
the photons as they are transported through the column via diffusion and advection, and 
ultimately escape through the column walls. We plot iV G as a function of the energy ratio 
e/eo and the location y in Figure 3 for the parameter values (3 = 0.4 and yo = 0.9. In this case 
the first eigenvalue is given by A = 4.231, and the corresponding photon index is a = 2.231 
(see Fig. 2). The selected value of y corresponds to a source located near the bottom of the 
accretion column, just above the stellar surface. At the source location, y = y = 0.9, the 
energy spectrum extends down to the injection energy, eo- However, at all other radii the 
spectrum displays a steep turnover above that energy because all of the photons at these 
locations have experienced Fermi energization due to collisions with the infalling electrons. 
The photons with energy e = eo at the source location have been injected so recently that 
they have not yet experienced significant energization. Note that for small values of y, the 
escaping spectrum is greatly attenuated due to the inability of the photons to diffuse up- 
stream through the rapidly infalling plasma. The escaping spectrum is also attenuated in 
the far downstream region (y ~ 1) due to the very high opacity of the flow, which inhibits 
the diffusion of photons through the walls. This latter effect stems from the divergence of the 
electron density near the stellar surface. The energy of maximum brightness, corresponding 
to the peak in iV G , achieves its maximum value as y — > because the photons that manage 
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Fig. 3. — Green's function N^(y ,y,e ,e) describing the photon distribution escaping from 
the accretion column per unit time (eq. [77]) plotted in units of A r ((T||/o-_|_) 1 / 2 (r e ) _1 as a 
function of the photon energy ratio e/eo for the indicated values of the spatial variable y. In 
this example we have set the absorption constant /3 — 0.4 and the source location yo = 0.9, 
so that the monoenergetic source is located near the base of the accretion column. 

to diffuse far upstream from the source are the ones that have resided in the flow the longest 
and therefore experienced the most energy amplification. 

We plot ATp as a function of energy and location in Figure 4 for the parameters /3 — 4 
and yo = 0.4, in which case we obtain for the first eigenvalue and the photon index Ao = 6.325 
and ao = 4.325, respectively. In this scenario the source is located in the upstream region 
and the absorption is stronger than that in Figure 3. Due to the increased absorption 
resulting from the larger value of f3, the photons spend less time on average in the flow being 
energized by collisions with the electrons before they escape through the column walls or 
are "recycled" by absorption. This causes a steepening of the spectrum at high energies, as 
evidenced by the increase in the photon index a . In this situation the energy of maximum 
brightness (where iV e G displays a peak) achieves its greatest value in the downstream region. 
This is the reverse of the behavior displayed in Figure 3 because in the present example, 
the source is located in the upstream region and therefore the photons that diffuse further 
upstream towards y = do not experience as much compression as those that are advected 
downstream. The escaping photon distribution in the far upstream and downstream regions 
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Fig. 4. — Same as Fig. 3, except (3 = 4.0 and y = 0.4. In this case the source is located 
in the upstream region, and the absorption at the source location is stronger than in Fig. 3. 
The latter effect causes a significant steepening of the spectrum at high energies, as explained 
in the text. 

is strongly attenuated due to the same processes operative in Figure 3, and therefore most 
of the escaping radiation is emitted from the accretion column around the source location. 

The analytical results for the Green's function obtained in this section provide the basis 
for the consideration of any source distribution since the fundamental differential equation 
(31) is linear. This is further discussed in § 5, where we use equations (35) and (44) to 
convolve the Green's function with the blackbody spectrum produced by the thermal mound. 



4.5. Column-Integrated Green's Function 

By integrating over the vertical structure of the accretion column, we can compute the 
total emitted radiation distribution, which corresponds to the phase-averaged spectrum of 
the X-ray pulsar. In the cylindrical geometry employed here, the formal integration domain 
is the region — oo < x < x st , where x st is the distance between the sonic point and the stellar 
surface (see eq. [21]). For the case of a monochromatic source, we define the column-integrated 



-25 - 



Green's function by writing 



/Est 
iV e G (a;o,x, e ,e) <ix , 
-oo 



(79) 



where $ G <ie represents the number of photons escaping from the column per unit time 
with energy between e and e + de. Using equation (50), the variable of integration can be 
transformed from x to y to obtain the alternative form 



<5>?(yo,e ,e) 



r 
2V3 



1/2 »i 



N?(y ,y,e ,e) 



dy 

y 



(80) 



where yo is related to x via equation (52). Despite the appearance of the factor y~ l inside the 
integral in equation (80), the contribution from small values of y is actually negligible because 
the spectrum declines exponentially in the upstream region due to advection. Furthermore, 
the escaping spectrum is also strongly attenuated in the downstream region due to the 
divergence of the electron density. Hence most of the radiation is emitted from the column 
around y ~ 0.5, as indicated in Figures 3 and 4. 

By combining equations (77), (78), and (80), we can reexpress the column-integrated 
Green's function as 



$f(yo,eo,e) 



N 

3/4 

e Uq n=0 



E 



9n(yo) 



3 — Xn 



(81) 



where the quadratic normalization integrals I n are evaluated using equation (72) and we 
have made the definition 



X n = g n (y) (l-y)y 1 dy . 
Jo 



(82) 



In Appendix E, we demonstrate that the integral in the expression for X n can be carried out 
analytically to obtain 



X n — 



+ 
+ 



T(l-a)(l-y ) 2 (yt c MK,y ) 



r(a) y? 2 (A„,y ) 

yo 



r(3 - c) 



3-c 

(c-2) 



F(3-a,3-6;4-c; y ) 
F(a 



+ 



F(2-a,2-6;3-c; y ) 



(a -2) (6-2) 



(a-l)(6-l)r(l-c 
'^- 2 ^ 0) ]} + (a-V-l) 



1 + 



F(a,6;c- l;y ) 

(c-2) 
(a -2) (6-2) 



(83) 



where v?i(A n ,y ), <fl(K,yo), and v? 2 (A„,y ) are given by equations (59), (60), and (62), 
respectively, and the constants a, b, and c are computed in terms of \ n using equations (61). 
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In Figure 5 we plot the dependence of the column-integrated Green's function <&f on 
the radiation energy e using the same parameters employed in Figures 3 and 4. Note that 
when (3 = 0.4 and yo — 0.9, which corresponds to the thin line in Figure 5, the column- 
integrated spectrum displays a peak at e/e ~ 6. The peak forms because the source is 
located close to the bottom of the accretion column, whereas most of the photons escape 
from higher altitudes around y ~ 0.5 (see Fig. 3). The escaping radiation has therefore 
experienced significant compression in the flow and consequently the photon energies are 
generally boosted well above the injection energy e . On the other hand, when (3 = A 
and yo = 0.4, which corresponds to the thick line in Figure 5, the absorption is so strong 
that the spectrum is quite steep and therefore the peak is suppressed. Consequently the 
column-integrated spectrum achieves it maximum value at the injection energy in this case. 



4.6. Reprocessed Blackbody Radiation 

The closed- form solutions for the Green's function iV e G (eq. [77]) and for the column- 
integrated Green's function <3>p (eq. [81]) provide a very efficient means for computing the 
X-ray spectrum escaping through the walls of the accretion column due to the injection of 
monochromatic seed photons. However, in our physical application to X-ray pulsars, we 
are primarily interested in computing the emitted spectrum resulting from the reprocessing 
of blackbody radiation injected into the accretion column from the thermal mound. In the 
escape probability approach utilized here, the photon distribution emitted through the walls 
of the cylindrical column at position x, denoted by N e (xo,x,e), is computed using 

2 2 

N e (x , x, e) = f(x , x, e) , (84) 

^esc 

where the particular solution f(x ,x,e) is evaluated using equation (44), and N t dxde rep- 
resents the number of photons escaping from the column per unit time between positions x 
and x + dx with energy between e and e + de. By combining equations (44) and (84), we 
find that the escaping photon spectrum can be expressed as 

N,(x , x, e) = ^ f fo(^e ,e) ^ ^ ^ (g5) 

£esc JO Jvo 

where the blackbody source function S'(eo) is evaluated using equation (35). The spatial 
variables (x,x ) and (y,yo) are interchangeable by virtue of equations (47) and (52), and 
therefore we can use equations (75) and (85) to write 

Ne(yo,y,e)= e S{e )de , (86) 

Jo 
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Fig. 5. — Column-integrated Green's function <3>^(y , e , e) describing the photon spectrum 
escaping through the walls of the accretion column (eq. [81]) plotted in units of iV /e as a 
function of the photon energy ratio e/eo- The parameter values used are yo = 0.9, (3 = 0.4 
(thin line) and yo = 0.4, (3 = 4.0 (thick line), which correspond to the spectra plotted in 
Figs. 3 and 4, respectively. 

where N^(y ,y,eo,e) is computed using equation (77). 

By analogy with equation (79) for the Green's function, we define the particular solution 
for the column-integrated photon spectrum escaping from the plasma due to the blackbody 
source, $ e (x ,e), by writing 



so that <3> e de gives the number of photons escaping from the column per unit time with 
energy between e and e + de. The variable of integration can be transformed from x to y 
using equation (50), which yields 



where y is related to x by equation (52). By substituting for N t using equation (86) and 




(87) 




(88) 
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employing equation (80), we can obtain the alternative form 

Hvo, e) = f ^ ( ^ ; e0 ' e) e 2 S(e ) de . (89) 
Jo iVo 

The column-integrated Green's function $f(y Q , €q, e) is evaluated using equation (81). In § 5 
we will utilize equations (86) and (89) to compute the photon spectrum emitted from the 
accretion column using parameters corresponding to specific X-ray pulsars. 

5. ASTROPHYSICAL APPLICATIONS 

For given values of the stellar mass M* and the stellar radius R*, our model has three free 
parameters, namely the column radius r , the temperature at the top of the thermal mound 
To, and the accretion rate M. In this section we investigate the relationship between these 
quantities and the dimensionless parameters yo and f3 appearing in the theory. Although it 
is not our intention in this paper to develop complete models for X-ray pulsar spectra, it is 
nonetheless interesting to compare the simplified model developed here with actual data for 
a few sources. Hence in this section we will also compute the spectrum emerging from the 
accretion column for parameters corresponding to two specific X-ray pulsars. 

5.1. Location of the Thermal Mound 

In order to understand how the theory developed here can be related to observables 
such as the temperature and the luminosity, we first need to determine how the position 
of the thermal mound, x Q (or, equivalently, y ), depends on the values of r , T , and M. 
The blackbody surface of the mound represents the photosphere for photon creation and 
destruction in this problem, and therefore the opacity is dominated by free-free absorption 
inside the mound. For a given value of the temperature To, the density at the top of the 
mound, p ? can be calculated by setting the free-free optical thickness of the column equal 
to unity, so that 

r ff = r a|(ro) = l, (90) 

where the Rosseland mean of the free-free absorption coefficient is evaluated in cgs units 
using (Rybicki & Lightman 1979) 

al(r ) = 6.10 x 10 22 T~ 7 ' 2 p 2 , (91) 
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for pure, fully-ionized hydrogen with the Gaunt factor set equal to unity. The density at the 
top of the thermal mound is therefore given by 

p = 4.05 x 1(T 12 T 7/4 r 1/2 . (92) 

The velocity at the top of the mound, v , is related to p via the continuity equation M = 
nrl v po, so that 

v = 7.86 x 10 10 Mr 3/2 T _7/4 . (93) 

By combining this result with equations (24) and (53), we find that the value of y is given 
by 

y = 1 - 2.15 x 10 14 R l J 2 M; 1 ' 2 M r" 3/2 T " 7/4 , (94) 

or, equivalently, 

The values obtained for y are extremely close to unity, indicating that the top of the 
thermal mound is just above the stellar surface. Once yo is computed, we can determine xq 
by using equation (52) to write 

^r 1+ wm' (96) 

where x st is the distance between the sonic point and the stellar surface given by equa- 
tion (21). The height of the thermal mound above the stellar surface is then given by (see 
eq. [22]) 

ho = h(x ) = x st , (97) 

or, equivalently, 

where we have substituted for x st using equation (21). Specific numerical examples are given 
in § 5.4. For typical X-ray pulsar parameters, we find that the mound is situated very close 
to the surface of the star, as expected. 



5.2. Computation of the Eigenvalues 



In order to compute the eigenvalues A n using equation (67), we must also evaluate the 
dimensionless absorption parameter f3 appearing in the transport equation (31). Based on 
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energy conservation considerations in the radiation-dominated plasma, we have shown in 
§ 3.2 that the photon energy flux Q must be continuous across the location of the thermal 
mound. This leads to equation (43), which can be rewritten as 

0=^' (") 

where P — U /3 is the radiation pressure at the top of the thermal mound. According to 
equation (49), P = (7/4) Jv c y , and therefore 

(3= AaT ° . (100) 
21 Jv c v i/o 

By substituting for v c and Vo using equations (24) and (53), respectively, and setting J = 
M/(7T7o), we can obtain the alternative form 

P= ■ , (101) 



6y (l-y )GM*M ' 



or, equivalently, 



2.24x10-3 / R. \(M*Y[ M (J^Y(Jh_X (102) 

P 2/o(l-l/o) VlOkm^ \M Q J ^lO^gs-y vlkm/ \WkJ ' 1 1 

Taken together, equations (95) and (102) allow the determination of the two dimensionless 
model parameters y and (3 in terms of R*, M*, r , T , and M. This closes the system and 
facilitates the calculation of the emergent photon number spectrum N e using equation (86), 
and the calculation of the column- integrated spectrum $ e using equation (89). The values 
of the parameters r , T , and M can therefore be extracted by comparing the theoretically 
predicted photon count-rate distributions with X-ray pulsar spectral data. 



5.3. Similarity Variables 

The eigenvalues \ n are functions of the dimensionless parameters y and (3. Based on the 
dependences of these parameters on r , T , and M expressed by equations (95) and (102), we 
note that the eigenspectrum remains unchanged if r , T , and M are simultaneously varied 
according to the prescription 

r oc T ' 9/2 oc M 9 / 10 . (103) 

In this case, the column-integrated Green's function (eq. [81]) also remains unchanged since 
it depends only on the values of f3 and y . These findings suggest that it is useful to introduce 
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the new "similarity variables" p and q, defined by 

^(iw)(io^) ' g " (wk) (irk) 2 ' 9 • (104) 

In terms of these parameters, equations (95) and (102) for y and (3 now become 

» = ^"^(vL) V2 {wj m ^ ^ • < 105 > 

and 

2.24 x 10- 3 [ R, \ I M.y 1 _ 5 9 

The introduction of the new variables p and q is convenient because it reduces the dimension- 
ality of the parameter space that determines the eigenvalues A n from the three-dimensional 
domain (r ,T ,M) to the two-dimensional space (p,q). In Figure 6 we plot the spectral 
index a — — 2 of the photon count-rate spectrum as a function of p and q, where Ao 
is the first eigenvalue and we have set R* = 10 km and M* = 1.4 M Q (the same quantity 
was plotted as a function of yo and (3 in Figure 2). Note that when q is increased for a 
fixed value of p, the values of yo and (3 both increase according to equations (105) and (106). 
The source therefore moves further downstream, and the absorption at the source location 
becomes stronger. This causes the spectrum to steepen as we move upwards along a vertical 
line in Figure 6. However, if q becomes very large, the spectrum starts to harden again due 
to the increasing amount of compression experienced as the source approaches the base of 
the accretion column. 



5.4. Example Spectra 

Our primary goal in this paper is to explore the effect of bulk Comptonization on the 
shape of the X-ray continuum spectrum emerging from a pulsar accretion column. The devel- 
opment of complete models that include additional effects such as thermal Comptonization, 
cyclotron processes, free-free emission and absorption, and iron line formation will be pur- 
sued in future work. Here we compare spectra computed using our simplified model with the 
observations for a few sources in order to focus attention on the role of the accretion shock 
in the formation of the X-ray continuum. In Figures 7 and 8 we plot the position-dependent 
photon count rate spectra N e (eq. [86]) escaping from the accretion column computed using 
the two sets of model parameters listed in Table 1. These results were obtained using 20 
eigenvalues and eigenfunctions, which yields extremely high precision due to the rapid con- 
vergence of the expansion. The theoretical spectra describe the upscattering of Planckian 
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Fig. 6. — Contour plot of the high-energy power-law photon spectral index cto as a function 
of the similarity variables p and q, which are related to the fundamental physical parameters 
r , T , and M via eqs. (104). The value of a is indicated for each curve, and we have set 
= 10 km and M* = 1.4 M Q . The parameters used for models 1 and 2 are represented by 
the labeled points. 

seed photons injected into the base of the accretion column from the surface of the thermal 
mound. In both cases we use M* = 1.4 M and i?* = 10 km for the stellar mass and radius, 
respectively. For model 1 (Fig. 7), the parameters adopted are r = 6 km, T = 7.3 x 10 6 K, 
and M = 2.69 x 10 16 gs _1 . In model 2 (Fig. 8) we set r = 1.3 km, T = 9.0 x 10 6 K, and 
M = 3.23 x 10 13 gs -1 . We use equations (104), (105), and (106) to compute the theoretical 
parameters p, q, yo and (3 based on the selected values for ro, Tq, and M. In model 1 we 
obtain (5 = 26.4505 and y = 0.99977, and in model 2 we find that (3 = 2.894 x 10 5 and 
y = 0.999998. Additional parameters are provided in Table 1. The fact that the thermal 
mound location y is extremely close to unity in both cases indicates that the mound lies just 
above the surface of the star, as expected. It is interesting to compute the actual physical dis- 
tance between the thermal mound and the stellar surface, ho, which is given by equation (98). 
Using the parameters corresponding to models 1 and 2 yields ho = 39.84 (ox/cry) 1 / 2 cm and 
h = 0.075 ((t±/ct\\) 1/2 cm, respectively. Note that in model 2 the thermal mound is essen- 
tially on the surface of the star, whereas in model 1 there is a much larger separation. This 
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Fig. 7. — Escaping photon distribution N e (yo, y, e) (eq. [86]) divided by (cr\\/(J±) 1 / 2 plotted as 
a function of e for the indicated values of y. These results were computed using the model 1 
parameters and describe the spectrum emitted per unit length of the accretion column due 
to the bulk Comptonization of blackbody radiation. 

reflects the fact that the accretion rate in model 1 is over 800 times larger than in model 2. 

The results obtained for the first eigenvalue and for the photon index in model 1 are 
Ao = 4.0398 and a = 2.0398, respectively. In model 2, we find that the spectrum is 
significantly steeper, with Ao = 4.6382 and ao = 2.6382. The steeper spectrum obtained in 
model 2 results from the large value obtained for /3, as indicated in Table 1. In general, the 
escaping spectrum has a power-law shape at high energies, as expected for a Fermi process, 
and a low-energy turnover due to the Planck distribution of the seed photons. Note that the 
emitted radiation is concentrated in a layer just above the thermal mound location, around 
y ~ 0.5. This reflects the fact that advection due to collisions with the infalling stream of 
high-speed electrons tends to keep the photons trapped at low altitudes in the column, and 
therefore few of them are able to escape through the column walls for y <C 1. Since the 
spin axis and the column axis are not aligned in pulsars, the portion of the accretion column 
visible along the line of sight to the Earth changes as the pulsar rotates. When combined 
with the vertical variation of the energy dependence of the escaping photon distribution, this 
effect will produce a pulse-phase dependence in the observed X-ray spectrum. 
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Fig. 8. — Same as Fig. 7, except the results were computed using the model 2 parameters. 
In Figures 9 and 10 we plot the photon count rate flux measured at Earth, 

Fe(e) = ~£kLP~ ' (1 ° 7) 

computed using models 1 and 2, respectively, where D is the distance to the source and 
the column-integrated spectrum $ e (y ,e) is evaluated using equation (89). Also included 
in Figures 9 and 10 for comparison are plots of the unfolded phase-averaged spectra for 
the X-ray pulsars 4U 1258-61 (GX 304-1) and 4U 0352+30 (X Persei), respectively. The 
4U 1258-61 data was reported by White et al. (1983), and the 4U 0352+30 data is the 
result of XSPEC analysis of an archival RXTE observation taken in July 1998 and reported 
by Delgado-Marti et al. (2001). Based on estimates from Negueruela (1998), the values 
adopted for the distance D are 2.5 kpc for 4U 1258-61 and 0.35 kpc for 4U 0352+30. The 
observational data represent the de-convolved (incident) spectra, which depend weakly on 
the detector response model. The theoretical spectra were computed using 20 eigenvalues 
and eigenfunctions for high accuracy, and various amounts of interstellar absorption were 
included as indicated in the figures. 

Although the results presented here are not fits to the data, we note that the general 
shape of the pulsar spectrum predicted by the theory agrees fairly well with the observations 
in each case, including both the turnover at low energies and the power law at higher energies. 
Several other sources yield similar agreement. In our model, the turnover at ~ 2keV is due 
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to Planckian photons that escape from the accretion column without experiencing many 
scatterings. This effect will tend to reduce the amount of absorption required to fit the 
observational data, compared with the amount required using the standard ad hoc models 
usually employed in X-ray pulsar spectral analysis (see the discussion in Hickox, Narayan, 
& Kallman 2004). 

The second source, X Persei, presents an interesting test for the model due to its rel- 
atively low luminosity, Lx ~ 10 34 ergs s -1 . The agreement between the theory and the 
observations suggests that radiation pressure may still be playing an important role in the 
dynamics of this source, despite its low luminosity. In fact, Langer & Rappaport (1982) 
pointed out that radiation pressure may have a strong effect on the flow dynamics of low- 
luminosity pulsars if the photon energies are close to the cyclotron resonance, which greatly 
increases the electron scattering cross section. Furthermore, the dynamical importance of 
the gas pressure in this luminosity range may lead to the development of a gas-mediated 
"subshock" within the overall gradual compression of the radiation-dominated shock. This 
type of hybrid shock is analogous to the "two-fluid" model of cosmic-ray mediated supernova 
shocks (e.g., Drury & Volk 1981). It is know that in the cosmic-ray case, multiple flow solu- 
tions can occur for the same set of upstream boundary conditions (e.g., Becker & Kazanas 
2001). There is some possibility that the same type of multi-state behavior may occur in 
X-ray pulsars, which could have interesting observational implications. 



6. MODEL SELF-CONSISTENCY 

The analytical approach employed in this paper is based on the conservation of energy 
and momentum in the radiation-dominated flow. This leads to the self-consistent determi- 
nation of the thermal mound location, xq, for given values of the fundamental parameters 
r , T , and M, as discussed in § 5.1. In this section we verify that the model satisfies the 
self-consistency conditions required for global energy conservation. We also examine the 
implications of the one- dimensional escape probability formalism that has been utilized to 
model the emission of radiation from the column. 
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Fig. 9. — Theoretical column-integrated count-rate spectrum F e (e) (eq. [107]) computed us- 
ing the model 1 parameters along with various amounts of interstellar absorption is compared 
with the X-ray spectrum of 4U 1258-61. The column densities used are Nn = (solid line); 
Nh — 3 x 10 21 cm~ 2 (dashed line); and JV H = 9x 10 21 cm~ 2 (dot-dashed line). The other 
parameters for the theoretical model are given in Table 1. 



The exact dynamical solution for the radiation pressure inside the pulsar accretion 
column is given by equation (29), which states that 



where x st is the distance between the sonic point and the stellar surface, given by equa- 
tion (21). In principle, the radiation pressure can also be computed by integrating the par- 
ticular solution for the photon distribution / inside the accretion column. Since P = U/3, 
we can use equation (32) to write 

1 f°° 

P(x) = - J e 3 f(x ,x,e)de, (109) 
where / is computed using the convolution (see eq. [44]) 



6.1. 



Pressure Distribution 




(108) 




(110) 
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Fig. 10. — Same as Fig. 9, except the data corresponds to 4U 0352+30 and the theoretical 
spectra were computed based on model 2. The column densities used are iV H = (solid 
line)] jV H = 3 x 10 21 cm" 2 (dashed line)] and jV H = 9x 10 21 cm" 2 (dot-dashed line). 

Here, the blackbody source S'(eo) is given by equation (35) and the Green's function f G is 
evaluated using equation (68), bearing in mind that the quantities (x,xq) and (y,yo) are 
interchangeable via equations (47) and (52). 

We can check the accuracy of the entire solution technique, including the determination 
of the eigenfunctions and the eigenvalues, by comparing the exact dynamical result for P 
with that obtained by integrating the spectrum. In Figure 11 we plot the exact analytical 
solution for the radiation pressure computed using equation (108) along with the results 
obtained by integrating the model 1 and 2 spectra using equation (109). Note that the 
agreement between the various solutions for P is excellent, and in each case the pressure 
achieves the correct stagnation value P st = (7/4) Jv c (eq. [30]) as the gas approaches the 
stellar surface. The agreement between the various results for the pressure confirms that our 
solution for the radiation spectrum is consistent with the dynamics of the accretion flow. 
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Fig. 11. — Radiation pressure P plotted in units of Jv c as a function of x, which is the 
distance measured from the sonic point in the direction of the flow. The dashed vertical line 
at x = represents the sonic point, and x = x st at the stellar surface. The exact dynamical 
solution (eq. [108]) is denoted by the solid line, and the filled circles and triangles represent 
the values obtained using eq. (109) to integrate the model 1 and 2 spectra, respectively. At 
the surface of the star, the pressure approaches the stagnation value (eq. [30]) as expected. 



6.2. Total Luminosity 

In order to further explore the self-consistency of the formalism developed here, we have 
also computed the total X-ray luminosity, Lx, emerging from the column by integrating the 
escaping number spectrum $ e (eq. [87]) using 

/•OO 

L x = e$ e (x ,e)de . (Ill) 
Jo 

We have confirmed that the values obtained using this expression agree precisely with the 
accretion luminosity given by 

GM*M , , 

Lx = • (112) 

The calculations performed in this section and in § 6.1 validate our entire solution procedure 
for the determination of the radiation spectrum by demonstrating that our results are fully 
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consistent with global energy conservation. 



6.3. Escape-Probability Formalism 

Within the context of our one-dimensional model, the diffusion of radiation through 
the walls of the accretion column is treated in an approximate manner using an escape- 
probability formalism. In this approach, the probability per unit time that a randomly- 
selected photon at location x inside the column will escape through the walls is equal to t~ s l, 
where the mean escape time t csc is a function of x given by (see eqs. [12]) 

tUx) = , (113 ) 

c 

with r and a± denoting the radius of the cylindrical column and the electron scattering 
cross section for photons propagating perpendicular to the column axis, respectively. The 
escape of radiation through the walls of the accretion column can therefore be modeled 
approximately using a term of the form 

df 
dt 



= -/-, (H4) 

approx esc 



which is incorporated into the transport equation (31). 

From a technical point of view, the utilization of the escape-probability formalism im- 
plies that the distribution of the photon number density across the column is proportional 
to the first eigenfunction of the operator describing the diffusion of radiation perpendicular 
to the column axis (e.g., Sunyaev & Titarchuk 1980). In a comprehensive, multidimensional 
calculation, the diffusion of photons perpendicular to the column axis is treated by replacing 
the escape-probability term df /dt = —f/t csc in the transport equation (31) with the more 
rigorous expression 



19 f cr df\ 



dt exact r Or \3n e a ± Or 
where r is the radial coordinate measured from the central axis of the cylindrical accretion 
column. If the gas density is constant across the column at a given height, as assumed here, 
then equation (115) can be separated to obtain for the first eigenfunction the solution 

/ ( r , t ) = x^,(-^) Jb (g, aw) 

where A is a constant, £ denotes the first eigenvalue, J represents the zeroth-order Bessel 
function, and 

e± = (n e (7 ± )- 1 (117) 
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is the electron scattering mean free path for photons propagating perpendicular to the column 
axis. 

The value of £ is determined by satisfying the flux boundary condition at the surface 
of the column. By employing the Eddington approximation, we can write the boundary 
condition as 

?±%+lf =0, (118) 

OT r=ro 

where 7 is a positive constant of order unity and we have assumed that no radiation is 
incident on the column from the outside. The precise value of 7 is somewhat arbitrary, 
with Rybicki & Lightman (1979) setting 7 = 3 1/2 and Sunyaev & Titarchuk (1980) setting 
7 = 3/2. Combining equations (116) and (118), we find that £ is the smallest positive root 
of the equation 

* m - © * m < 

where J\ denotes the order-unity Bessel function. 

The rate of change of the first eigenfunction (eq. [116]) is given by 

of c e 



dt ' 3£i 



f(r,t). (120) 



Comparison of equations (114) and (120) reveals that the escape-probability approximation 
and the rigorous diffusion model yield the same result for the logarithmic time derivative 
dhif/dt if the mean escape time t eS c and the first eigenvalue £ are related via 

tesc = ^ • (121) 

This equation is well satisfied in all of our sample calculations, which suggests that the 
escape-probability formalism provides a reasonable description of the leakage of photons 
from the accretion column. 



7. CONCLUSIONS 

We have shown that the simplified model considered here, comprising a radiation- 
dominated accretion column with a blackbody source/sink at its base and a radiative shock, 
is able to reproduce the power-law continuum spectra observed in accretion-powered X-ray 
pulsars with a range of luminosities, as indicated in Figures 9 and 10. Our results represent 
the first ab initio calculations of the X-ray spectrum associated with the physical accretion 
scenario first suggested by Davidson (1973). For given values of the stellar mass M* and 
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the stellar radius R*, our model has only three free parameters, namely the column radius 
r , the temperature at the top of the thermal mound T , and the accretion rate M. We 
have confirmed that the solution obtained for the radiation spectrum is consistent with the 
dynamics of the flow in the accretion column (see Fig. 11). Many pulsars are found to have 
spectra that are well fitted using a combination of a power-law plus a blackbody compo- 
nent, but until now this form has been adopted in a purely ad hoc manner. The new model 
described here finally provides a firm theoretical foundation for this empirical result. Our 
analytical solution, based on a rigorous eigenfunction expansion, is relatively straightforward 
to implement computationally because it avoids numerical integration of the transport equa- 
tion. The work presented here therefore represents a significant step in the development of 
a comprehensive theory for the spectral formation process in X-ray pulsars. 

The photon spectral index a is related to the first eigenvalue A via a = A — 2, and 
therefore the leading eigenvalue determines the slope of the high-energy portion of the X-ray 
spectrum. The value of «o is quite sensitive to the strength of the absorption at the thermal 
mound, which is determined by the parameter f3. According to equation (102), (3 in turn is 
a strong function of the mound temperature To. The higher the temperature, the stronger 
the emission and absorption at the blackbody surface of the mound. Model 2 has a higher 
temperature than model 1, and therefore the value of f3 is larger. Hence the spectrum is 
steeper in model 2 because the photons spend less time on average upscattering in the flow 
before either escaping through the column walls or being "recycled" by absorption at the 
mound location. Acceptable values for a , in the range between 2 and 4, are obtained for 
reasonable source parameters, as indicated in Figures 2 and 6. Although a broad range of 
spectral indices can be produced, we must have a > 2 in order to avoid an infinite photon 
energy density since the model considered here does not include a high-energy cutoff. 

The shape of the Green's function spectrum derived here is qualitatively similar to 
the results obtained using the "bulk motion Comptonization" model coded into the XSPEC 
package using the BMC function (see Shrader & Titarchuk 1998, 1999; Borozdin et al. 1999). 
The BMC model includes parameters corresponding to the high-energy spectral index, the 
temperature of the soft photon source, the fraction of the Comptonized component in the 
resulting spectrum, and the size of the emitting region. However, the BMC model is designed 
to treat spectral formation in black-hole accretion flows, whereas our focus here is on neutron 
star accretion. Hence the model developed in this paper can be viewed as a parallel formalism 
that is optimized for the treatment of the spectral formation process in accretion-powered 
X-ray pulsars. In this sense, our theory translates the free parameters of a generalized bulk 
Comptonization model into specific values for the accretion rate, temperature, and column 
radius for a particular X-ray pulsar. 
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Our goal in this paper is to explore the direct role of the accretion shock in producing the 
observed continuum spectra in X-ray pulsars via the first-order Fermi process. Although the 
model produces spectra that are quite similar to those observed from several X-ray pulsars, it 
cannot explain the spectra from other sources. For example, the observed spectra for several 
of the brightest pulsars with luminosities in the range 10 37 ™ 38 ergss" 1 (e.g., Her X-l) are very 
flat power-laws leading up to a sharp exponential cutoff. The photon index in these sources 
is typically a ~ 1, which is outside the range allowed by our model since it does not include 
a high-energy cutoff. In order to understand spectral formation in the brightest pulsars, 
additional effects such as thermal Comptonization, cyclotron features, bremsstrahlung, and 
iron emission must be incorporated into the model. In principle, one should also consider 
higher-order relativistic effects that can alter the form of the fundamental transport equation 
(Psaltis & Lamb 1997). Thermal Comptonization is particularly important since this pro- 
cess redistributes energy from the highest frequency photons to lower energies via electron 
recoil. This naturally leads to the development of a high-energy exponential turnover and a 
concomitant flattening of the spectrum at moderate energies, as observed. 

The thermal Comptonization process was studied by Lyubarskii & Sunyaev (1982) and 
Becker (1988) in the context of plane-parallel shocks. However, the velocity profiles em- 
ployed by these authors are not consistent with the structure of an X-ray pulsar accretion 
column. Poutanen & Gierlihski (2003) computed X-ray pulsar spectra based on the thermal 
Comptonization of soft radiation in a hot layer above the magnetic pole, but their model 
did not include a complete treatment of the bulk process, which is of crucial importance in 
X-ray pulsars. In this paper we have carefully analyzed the effect of bulk Comptonization 
on the emitted spectrum, while neglecting the corresponding thermal process. The inclu- 
sion of thermal Comptonization is beyond the scope of the present paper because it renders 
the transport equation insoluble via analytical means. However, in future work we intend 
to develop a comprehensive numerical model that incorporates both the bulk and thermal 
processes. 

The authors would like to thank Drs. Lev Titarchuk, Philippe Laurent, and Nikos 
Kylafis who provided useful comments on the manuscript. The authors are also grateful to 
the anonymous referee whose careful reading and stimulating suggestions led to significant 
improvements in the paper. PAB would also like to acknowledge generous support provided 
by the Office of Naval Research. 



-43- 



APPENDICES 

In the following sections we provide further details regarding the various mathematical 
methods that form the foundation of our model. 



A. Fundamental Solution in the Downstream Region 

When y ^ y , the spatial separation functions g(X,y) satisfy the homogeneous differen- 
tial equation (cf. eq. [55]) 

^ (1 -^H~.)a; + I^H 9=0 ' (A1) 

which has the fundamental solutions 

ipi(\,y) = yF(a, b; c; y) , (A2) 

<pl(X,y) = y- 1/A F(a - 5/4, b - 5/4; 2 - c; y) , (A3) 

where F(a,b; c; z) denotes the hypergeometric function (Abramowitz & Stegun 1970), and 
we have made the definitions 

9 - V17+ 16 A , 9 + V17 + 16 A 9 /4J . 

a = o ) o = o 5 c == — . (A4) 

8 8 4 v y 

As explained in the discussion following equation (61), in the downstream limit the radiation 
pressure approaches the finite stagnation value given by equation (30). The spatial separation 
functions g(\,y) should mimic this behavior, and therefore we require that g approaches a 
constant in the limit y — > 1. This does not happen in general, but only if the separation 
constant A is equal to one of the eigenvalues A n , in which case we obtain the eigenfunction 

g n (y) = g(K,y) ■ (A5) 

The hypergeometric functions appearing in equations (A2) and (A3) can be evaluated at 
y — 1 using equation (15.1.20) from Abramowitz & Stegun (1970), which gives for general 
values of a, b, and c 

, 1X T(c)T(c-a-b) ,.„, 
F M ;e;l) = r L4__j . (A6) 

However, for the values of a, 6, and c in the current application, we find that 

c - a - b = , (A7) 
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and therefore the hypergeometric functions F(a, b; c; y) and F(a — 5/4, b — 5/4; 2 — c; y) 
each diverge in the downstream limit y — > 1. 

Based on the asymptotic behaviors of the hypergeometric functions appearing in the 
expressions for ipi and <£>*, we conclude that in the downstream region (y > y ), the eigen- 
function g n must be represented by a suitable linear combination of </?i and that remains 
finite as y — > 1. In order to make further progress, we need to employ equation (15.3.10) 
from Abramowitz & Stegun (1970), which yields for general a, b, and y 



-M/(6 + n)-ln(l-y)l(l-y)", (A8) 



where 

By making use of this expression, we note that the logarithmic divergences of the two func- 
tions ipi and ifil in the limit y — > 1 can be balanced by creating the new function 

s FwTTT^ - ^ • < A10 > 

which remains finite as y — > 1. Hence represents the fundamental solution for g n in the 
region downstream from the thermal mound located at y — yo- We can use the asymptotic 
behaviors of ipi and </?* as y — > 1 to conclude that 

7r[cot(7ra) + cot(7r6)] 

h ^ (A '^ = r(a)r(i-6) • (A11) 

Since the solutions </?i and </?2 are applicable in the upstream and downstream regions, 
respectively, the global expression for the eigenfunction g n is therefore 

S»M = i£ l(A ";? \ (A12) 
where the constant _B„ is evaluated using the continuity condition 

It follows from equations (All), (A12), and (A13) that the downstream value of g n is given 

by 

/ \ 7r[cot(7ra) +cot(7r6)] y?i(A„,?/o) 

lim g n (y) = 77 tt r , A14 

w-i r(a)r(l-6) ^ 2 (A n ,y ) 

where a and 6 are evaluated using equations (61). 
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B. Analytical Expression for the Wronskian 



The derivation of the eigenvalue equation (67) in § 4.2 makes use of the Wronskian of 
the two functions ipi and f 2 , defined by 



W{\,y) = if! — cp 2 



dy 



dy 



(Bl) 



It is useful to derive an analytical expression for W based on the differential equation (55) 
governing the two functions ipi and ip 2 , which can be rewritten in the self-adjoint form 



d_ 

dy 



i/4 /-, \ d(fi 
V {1 ~ y) dy-. 



+ 



A 



4y3/4 



y?-Ty? = , 



where 



T= 1-y 3(3v 5(y-y ) 



(B2) 



(B3) 



4y 7 /4 7v c y 3 / 4 

By applying equation (B2) to the function ip 2 and multiplying the result by ipi, and then 
subtracting from this the same expression with ipi and ip 2 interchanged, we obtain 



d_ 

dy 



y 1/4 (i - y) 



d(f 2 
dy 



(f2 



d_ 

dy 



,V4 



y) 



d(p i 
dy 



= 



which can be rewritten as 



sl/4(1 - a) |r + w '^[ ! ' 1/4 < 1 - ! '>] = °. 



where we have made use of the result 

dW 



dy 



d 2 ip 2 d 2 ip\ 



dy 2 



dy 2 



(B4) 



(B5) 



(B6) 



Equation (B5) can rearranged in the form 

dlnW d f 1/4 , 

which can be integrated to obtain the exact solution 

D(X) 



W{\,y) 



yV4 ( i _ y) ' 



(B7) 



(B8) 



where -D(A) is an integration constant that may depend on A but not on y. The exact 
dependence of D on A can be derived by analyzing the behaviors of the functions ipi and ip 2 
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in the limit y — > 0. For small values of y, we have the asymptotic expressions (Abramowitz 
& Stegun 1970) 

<pi -> y , y -»■ o , 

- -rR^)^ 1 ' 4 - »-». < B9 > 

We therefore find that asymptotically, 

W _> 5 r (^- fl ) y -V4 o . (BIO) 

4 r(a) r(2 - c) y ' y K ' 

Comparing this result with equation (B8), we conclude that 

D(X) = 5 - ra ~ a \ , (Bll) 
V ; 4 r(a)r(2-c) V ; 

and therefore the exact solution for the Wronskian is given by 

This result is used in § 4.2 in the derivation of the eigenvalue equation (67). 



C. Orthogonality of the Eigenfunctions 



The calculation of the expansion coefficients C n in the series representation for the 
Green's function 



/g(J/o, V, eo, e) = e 3 Cn ( — 
„=o V e ° 



3 A-n 



(CI) 



discussed in § 4.3 depends on the establishment of the orthogonality of the eigenfunctions 
g n (y)- Here we shall demonstrate that the eigenfunctions form an orthogonal set as required. 
Let us suppose that g n (y) and g m (y) are two eigenfunctions corresponding to the distinct 
eigenvalues A n and A m , respectively. The functions g n and g m satisfy the differential equation 
equation (55), and therefore we can utilize the self- adjoint form to write (cf. eq. [B2]) 



d 



y 1/4 (i - y) ^ 

dy i 



+ 



4y 3 / 4 



g n - Tg n > = , 



and 



9n 



dy 



y l "{i-y) 



dg m 
dy 



+ 



4y 3/4 



9m ~ Tg r , 



, 



(C2) 



(C3) 
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where T is given by equation (B3). Subtracting the second equation from the first yields, 
after integrating by parts with respect to y from y — to y — 1, 

(An - A m ) f y^ /A g n {y) g m {y) dy = 4 y 1 " (1 - y) 
Jo 

Based on the asymptotic behaviors of the eigenf unctions g n and g m given by equations (A14) 
and (B9), we find that the right-hand side of equation (C4) vanishes exactly, and we therefore 
obtain 

(A„ - A m ) / y-^ 4 g n (y) g m (y) dy = , (C5) 
Jo 

which establishes the orthogonality of the eigenfunctions. 



9n 



dg m 
dy 



- gn 



dg n 
dy 



(C4) 



D. Quadratic Normalization Integrals 

The evaluation of the Green's function / G using equation (CI) relies on knowledge of 
the expansion coefficients C n , which are given by 



12 N y 3/i g n (y 
lnrlv c I n 

where the quadratic normalization integrals I n are defined by 



c n = -: oy \ , (di) 



in= f y- 3/A 9 2 n (y) dy . (D2) 
Jo 

The direct computation of the normalization integrals I n via numerical integration is costly 
and time consuming, and therefore it is desirable to have an alternative procedure avail- 
able for their evaluation. In fact, it is possible to derive an analytical expression for the 
normalization integrals based on manipulation of the fundamental differential equation (55) 
governing the eigenfunctions g n {y)- 

Let us suppose that g(X,y) is a general solution to equation (55) for an arbitrary value 
of A (i.e., not necessarily an eigenvalue) that has the asymptotic upstream behavior 

g(X,y)^y, y^O. (D3) 

We shall also require that g be continuous at y — y , and that it satisfy the derivative jump 
condition (see eq. [57]) 



,. dg 
lim — 

e^o dy 



y=yo+e 



dg_ 

dy 



y=yo-e 



^-9{\Vo). (D4) 
4 yo 
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After a bit of algebra, we find that the global solution for g can be expressed as 

(\ ) = ! V^'V) ' V<Vo , 

9{ > y> \(l + a) l p 1 (X,y) + b l p 2 (X,y) , y > y , 

where a and b are given by 

. = _ 3/fyi(A,s/o) (p 2 (\yo) - = 3/^(A,g/ ) 

4j/ W(A,y ) ' 4y W(A,yo) ' 

and the Wronskian W is evaluated using equation (B12). 



(D5) 



(D6) 



Comparing the general solution for g(X,y) with the solution for the eigenfunction g n (y) 
given by equation (A12), we note that 

lim a = -1 , lim b = B n . (D7) 

A — >\n A >\n 

We can now use the self-adjoint form of equation (55) to write (cf. eqs. [C2] and [C3]) 



9n 



d 
dy 



y 1/4 (i-y) 



dg_ 

dy 



+ 



A 



and 



9 



dy 



y l,A (i - y) 



dg n 
dy 



+ 



4 y 3/4 

A 



9-Tg 



4^3/4 



9n - T g n } = , 



(D8) 



(D9) 



Subtracting the second equation from the first and integrating by parts from y = to y = 1 
yields 



(A - A n ) f y-*" g(X, y) g n {y) dy = 4 y 1 " (1 - y) 
Jo 



i \ ,dg n dg 

9^y)^-9n{y) Wy 



(D10) 



Since g — > y and g n — > y as y — > 0, we conclude that the evaluation at the lower bound y = 
on the right-hand side yields zero, and consequently in the limit A — > A n we obtain for the 
quadratic normalization integral I n (see eq. [D2]) 



In= [ 1 y-Wfi(y)dy= lim 4yV4(1 ~ y) [9{Ky) ( d 9n/ d y)-9n(y) (dg/dy)] 

J0 A^A n A — X n 



[Dll) 



y=l 



The numerator and denominator on the right-hand side of equation (Dll) each vanish 
in the limit A — > A n , and therefore we can employ L'Hopital's rule to show that (e.g., Becker 
1997) 

'dg dg n d 2 g 



4= lim 4 2/ 1 / 4 (l-y) 

A^A n 



dy dy ® n dy dX 



(D12) 



y=l 



-49- 



Substituting the analytical forms for g n (y) and g(X,y) given by equations (A12) and (D5), 
respectively, we find that equation (D12) can be rewritten as 



lim Ay l l A (l-y)B n 



(D13) 



A— A^ 



where we have also utilized equations (Bl) and (D6). Based on the asymptotic behavior of 
</?2 given by equation (All), we conclude that the final two terms on the right-hand side of 
equation (D13) contribute nothing in the limit y — > 1, and therefore our expression for I n 
reduces to 

' lim Ay^(l-y)B n W(\,y) da 



dX 



(D14) 



A — X n 



Since y — 1 is a singular point of the differential equation (55), it is convenient to employ 
the relation (see eq. [B8]) 

W(X, y) y 1 ^ (1 - y) = W(X, y ) yT (1 - Vo) , (D15) 

which allows us to transform the evaluation in equation (D14) from y = 1 to y = yo to obtain 
the equivalent result 

^i(K,yo) dhia 



4 2 /o V4 (l 



y ) aW(X,y ) 



(D16) 



A — A n 



^2(X n ,y ) dX 

where we have also substituted for B n using equation (A13). The derivative on the right-hand 
side can be evaluated using equation (D6), which yields 

din a <91no3i dhupn <91nW^ 



dX = dX + ~dX dX~ ' 

where the derivative of the Wronskian is given by (see eqs. [A4] and [B12]) 

dlnW *(a) + *(l-a) 



dX 



and 



17+16A)V2 
1 dT(z) 



(D17) 



(D18) 



(D19) 



T(z) dz 

Combining equations (D6), (D16), (D17), and (D18), we find that that the quadratic nor- 
malization integrals can be evaluated using the closed-form expression 

I n = K(X n ,y ) , (D20) 

where 

'W(a) + W(l - a) dlnipi d\n(p 2 



K(X,y)=3f3y- 3 / 4 (l-y) V 2 1 (X,y) 



(D21) 



(17+16A)V2 dX dX 

This formula provides an extremely efficient alternative to numerical integration for the 
computation of I n . 
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E. Column-Integrated Green's Function 

In § 4.5 we demonstrated that the column-integrated Green's function $f (y ,eo,e) can 
be written as (see eq. [81]) 



eyJ „_ n ^ \eo/ 



y n \.yo) i i 3 " 

where the quadratic normalization integrals I n are computed using equation (72), and 



X n = / g„(y) (l-y)y 1 dy . 
Jo 



(E2) 



As an alternative to numerical integration for the evaluation of the integral X n , here we shall 
derive an analytical expression for this quantity. First we use equation (63) to substitute for 
g n (y) in equation (E2) to obtain 

x. = f"' MK, ») (i - ») * + ^4 /' MK, V) (1 - V) * . (E3) 
Jo V </?2(A„,y J Jy y 

Substitution for ip 2 using equation (62) and rearranging terms, we can write 

¥>i(A„,j/ ) f 1 i \ Wl vdj/ r ^!(A n ,y ) P° Wl ,dy 



X n = Li—- / Vl {X n ,y){l-y) L 2^\ r / <Pi (K, V) (1 - V) 

<P2{K,yo) Jo y <P2{K,yo) Jo 



y 



T <fi(X n ,y ) f 1 

— ^2 



where 



/ ^(A n ,y)(l-y)^, (E4) 

Jy ( > y 



Ll = T(c)T(l-b) ' L2 = r(2 - c) T(a) ' (E5) 
To make further progress, we need to evaluate the fundamental indefinite integrals 

K,= f cp 1 (\ n ,y)(l-y)^- , K 2 = I tf(A n , y) (1 - y) ^ . (E6) 

J y J y 

Fortunately, these integrals can be carried out analytically. By applying equation (15.2.1) 
from Abramowitz & Stegun (1970) twice, we obtain for K\ the result 

* = („ 1) (a - 2H6 - 2) - 2, » - % c - 2; v) 

+ (I-i> ^4y - F(a ~ 1 ' t ~ 1;c ~ 1;i,)+Cl ' (E7) 
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where C\ is an integration constant. Likewise, we can apply equation (15.2.4) from Abramowitz 
& Stegun (1970) twice to find that 



+ 



[i -y) y 

2-c 

y 3 ~ c 



2-c 



(2-c) (3-c) 



F(l-a,l-fe;3-c;y) 

F(l-a,l-6;4-c;y)+C 2 , 



(E8) 



where C 2 is another integration constant. Combining relations and simplifying, after some 
algebra we obtain the final result 



X n — 



+ 
+ 



r(l-a)(l-y ) 2 f y 2 - c <pi (Xn, Vo) 
r(a) <£ 2 (A n ,?/o) 

2/0 



r(3 - c ; 

F(3-a,3-6;4-c;j/ )l +7- 

J a 



F(2-a,2-6;3-c; y ) 

¥>l(Af»,J/o) 



3-c 
(a -2) (6-2) 



(a - 1) (6- 1)T(1 -c 
F ( a, & ;c-2; yo) ]} + (Q J 1 i) -^_ i) 



1 + 



F(a,b;c-l;y ) 

(c-2) 
(a -2) (6-2) 



(E9) 



This formula allows the efficient computation of the integrals X n appearing in equation (81) 
for the column-integrated Green's function. 
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Table 1. Model Parameters 



Model 


r (km) 


To (K) 


M (gs- 1 ) 


P 


1 


yo 








vo/c 


1 


6.0 


7.3 x 10 6 


2.69 x 10 16 


0.890 


1.087 


0.999770 


2.645 x 10 1 


2.0398 


2.71 x 10~ 4 


1.48 x 10~ 4 


2 


1.3 


9.0 x 10 6 


3.23 x 10 13 


0.286 


0.954 


0.999998 


2.894 x 10 5 


2.6382 


2.36 x 10" 6 


1.22 x 10~ 6 



